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Dynamic Nuclear Polarization (DNP) is to date the most effective technique to increase the 
nuclear polarization opening disruptive perspectives for medical applications. In a DNP setting, the 
interacting spin system is quasi-isolated and brought out-of-equilibrium by microwave irradiation. 

Here we show that the resulting stationary state strongly depends on the ergodicity properties of the 
spin many-body eigenstates. In particular the dipolar interactions compete with the disorder induced 
by local magnetic fields resulting in two distinct dynamical phases: while for weak interaction, 
only a small enhancement of polarization is observed, for strong interactions the spins collectively 
equilibrate to an extremely low effective temperature that boosts DNP efficiency. We argue that 
these two phases are intimately related to the problem of thermalization in closed quantum systems 
where a many-body localization transition can occur varying the strength of the interactions. 


Introduction. — Formulating statistical mechanics 
for isolated many-body system requires the tacit assump¬ 
tion of ergodicity. Only recently, however, the eigen¬ 
state thermalization hypothesis (ETH) [lj [2] promoted 
this concept to a testable condition at the quantum level, 
allowing the identification of situations where ergodicity 
might even be broken once disorder combines with quan¬ 
tum interference m- It is natural to ask, then, whether 
ETH influences also the stationary regimes where energy 
is constantly injected and dissipated, leading again to an 
emergent simple description. Dynamic Nuclear Polariza¬ 
tion (DNP), the most effective technique to increase the 
nuclear polarization, is a paradigmatic out-of-equilibrium 
protocol to test these ideas. In a DNP procedure [5], the 
compound is doped with radicals (i.e., molecules with 
unpaired electrons), exposed to a strong magnetic field 
at low temperature, /3 _1 , and then irradiated with mi¬ 
crowaves (see Fig. [l] for details). At thermal equilib¬ 
rium, the unpaired electrons are much more polarized 
than nuclear spins because the electron Zeeman gap is 
orders of magnitude larger than the nuclear one. When 
the microwaves are on, at a frequency close to the elec¬ 
tron Zeeman gap, the spin system of interacting electrons 
and nuclei organizes in an out-of-equilibrium steady state 
with a huge nuclear polarization. The hyperpolarized 
sample can then be dissolved at room temperature [6], 
injected in patients, and used as metabolic tracer [7]. 
However, our understanding of the physical mechanisms 
that trigger hyperpolarization is still poor. A striking 
experimental evidence is the thermal mixing of the en¬ 
semble of different nuclear spins ( 13 C, 15 iV, 89 F, ...) 
0 0: their enhanced polarizations are well described 
by an equilibrium-like polarization, P n = tanh(/3 s /kj n /2) 
(see Fig. [l] right). While the Zeeman gap cj n depends 
on the nuclear species, the spin temperature Pj 1 is a 
unique parameter, possibly one thousand times smaller 
than /3 _1 , the one of the bath [TO] . 

But how can a quantum system appear thermal and 
colder when irradiated by microwaves? In which way the 
spin temperature can be controlled acting on the exper¬ 
imental parameters? 



FIG. 1. Color online. A solid material containing nuclear 
spins (e.g. 13 C, 15 N) and doped with molecules with unpaired 
electrons (left). At 1.2 Kelvin and 3.35 Tesla the equilibrium 
polarization of the electron spins is very high, (94%), while 
nuclear spins are very little polarized, less than 1%. Under 
microwave irradiation the spin system evolves towards a new 
steady state characterized by a single spin temperature P~ x ~ 
1 mK (right). In this work, we analyze exclusively the electron 
spins and show that an out-of-equilibrium spin temperature 
results from the interplay of disorder and interaction. 


In this paper we show that the spin temperature con¬ 
cept is directly connected to quantum ergodicity and 
ETH. While for classical physics, thermalization has its 
origin at the onset of chaotic dynamics, quantum ergod¬ 
icity requires the ETH, a thermal behavior at the level 
of single eigenfunctions |lj [TT] . The realm of ETH is 
normally restricted to quench protocols in cold atoms 
experiments, where any exchange of energy with the en- 
viroment is under control. Our work shows that ETH 
may impose a thermal behavior to the stationary state 
of open quantum systems, giving a practical and exper¬ 
imental relevance to the fundamental problem of Quan¬ 
tum Thermalization 00]. 

The microscopic model. — The traditional descrip¬ 
tion of DNP in the thermal mixing regime relies on 
the phenomenological assumption that the electron spins 
cool down once irradiated and act as a reservoir for all 
nuclear species. Here, we focus only on the electron spins 
and on the origin of the spin temperature in their station- 
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ary state. In the electron spin Hamintonian, the presence 
of g-factor anisotropy induces a spread of the electron 
Zeeman gap: 


N 

H = hY / ^e + A,)S i z +H Aip . ( 1 ) 

i=1 


where N is the number of electrons and uj e is the exter¬ 
nal magnetic field. The random fields A i are quenched 
at low temperature and distributed according to the den¬ 
sity /(A), with mean A^ = 0 and variance A? = Ac< j 2 . 
The term Hdip contains the interactions between spins 
due to the dipolar coupling. Experiments can access 
the product f(A)P e (u: ), dubbed EPR spectrum, where 
P e ( uo) is the polarization of an electron with Zeeman 
gap uj = uj e + A^. At equilibrium with the environ¬ 
ment (sketched in blue in Fig. §, P e (w) ~ -1. The 
microwave irradiations at frequency cjmw and intensity 
takes the form H MW = 2ojiJ2iS x cos(o;mw£), with S x 
total spin operator along the x component. In absence of 
dipolar interactions, the Bloch equations predict that the 
electrons with a resonating Zeeman gap are saturated, 
P e (c u rsj ujmw) ~ 0, while the others remain highly po¬ 
larized. This corresponds to the hole burning shape of 
EPR spectrum, showed in Fig. [2jA. On the contrary, ac¬ 
cording to the thermal mixing picture, the presence of 
dipolar interactions induces a collective reorganization 
of the electron polarization profile P e (cj), that shows an 
equilibrium-like shape even under microwave irradiation 


P e (cj) = — tanh 


2 


(uj - UJ o) 


( 2 ) 


with cj 0 — The ansatz of Eq. Q lacks a micro¬ 

scopic derivation. Moreover, recent ab-initio models [12] 
have only observed an hole burning shape, with a weak 
polarization enhancement triggered by local hybridiza¬ 
tions [Mm. Here, we take 


^di P = E^(^-+ c - c -) • ( 3 ) 

i<j 


where the Aij are the dipolar couplings. Because of the 
glassiness of DNP samples, the distance between elec¬ 
trons is random and, thus, for simplicity the coupling 
are taken, within a mean field approximation, as gaussian 
random variables with zero mean and variance U 2 /N. 
We are interested in the strongly correlated regime where 
disorder and interaction compete, i.e. U ~ Auo e . Our 
conclusions should not depend on the specific model and, 
here, we choose a uniform distribution of local magnetic 
fields by taking equally spaced A i = Auo e ( 2z ~^~ 1 ), with 
randomness only affecting the dipolar couplings. 

The master equation. — A key experimental obser¬ 
vation mm is that the spin system is quasi-isolated 
with a dephasing time, X 2 , very short compared to the 



FIG. 2. EPR spectrum Equilibrium (blue) versus MW ir¬ 
radiated (yellow) profile. Under irradiation two possible pro¬ 
files are expected: (A) the hole burning shape, characteristic 
of the non interacting case (B) the hyperbolic tangent shape 
characterized by a very low effective temperature, /3 S , that 
cools down the nuclear spins. 


time-scales of microwave dynamics and to the relaxation 
time, Ti, with the thermal bath m • Therefore, any ini¬ 
tial density matrix, p is quickly reduced by dephasing to a 
diagonal form in the basis of eigenstates of H. In practice 
m, the Lindblad equation p = Cp used to describe the 
dynamics of the open system reduces to a master equa¬ 
tion for the time evolution of the 2 N occupation probabil- 
ities, p nn with rates W n ^ n , = h(Ae n , n ,)W*%Y + W™™, 
where 


<n' h = |-E E \(n\Si\n')\ 2 , (4a) 

1 j=l a=x,y,z 

ti^mw _ 4cj 2 T 2 | (n\ S x | n') | 2 

n ’ n ' 1 + T 2 (| Ae nn /1 — cj M w) 2 1 j 


Here the index n label eigenstate of energy e n with 
Ae nn / = e n — e n f. The function h(x) = e^ x /(l + e^ x ) 
assures the convergence to Gibbs equilibrium when the 
microwaves are off and the rate in (Eal comes from 


single spin flip transitions on a time scale T\. Eq. (4b) 
describes transitions induced by the microwave field. In 
Fig. [3] we present the stationary value of the polarization 
P e (cj = uo e + A^) = 2Tr[S'*Poo], computed from the sta¬ 
tionary occupation probabilities which solve £poo = 0. 
Note that this requires the full diagonalization of P, 
strongly constraining the possible system sizes. 


Two possible behaviors are observed: For weak inter- 












3 



FIG. 3. Electron polarization under MW irradiation with 
wmw — 93.8775 Ghz, for the model of Eq.|l]) with N = 12 
spins. For strong dipolar interactions (circles and diamonds) 
the spin temperature is well defined. The fit according to Eq. 

(red and green lines) gives f3 s 1 =3.7 mK for U = 15 Mhz 
and /3- 1 = 7.4 mK for U — 45 Mhz. For weak interactions 
U — 2 Mhz (square) a simple broadening of the hole burning 
non interacting profile (dashed blue line) given in Eq. (8) of 

m- 

actions, the hole burning shape, already observed in [H- 
□3, is recovered. Instead, in presence of strong dipolar 
interactions, we show that all electrons rearrange accord¬ 
ing to the spin temperature profile of ©• Remarkably, 
in the wing closer to cjmw? electron polarization can even 
invert its sign. 

The origin of these two dynamical regimes can be un¬ 
derstood in relation with the quantum thermalization 
of the electron spins. In general, in closed quantum 
systems, an arbitrary initial state converges to a time- 
independent density matrix because of dephasing: in the 
basis of eigenstates, off-diagonal elements are suppressed 
while the occupation probabilities on the diagonal remain 
constant. But then, how could thermodynamics emerge 
if the initial occupancies are conserved? In Fig. [4] (left) 
we report the most probable value of the local polariza¬ 
tion for the eigenstates at a given energy. In presence of 
weak interaction (U = 2 Mhz) the polarization fluctuates 
between the extremal values ±1 showing that the exact 
eigenstates are almost factorized on local spins. When 
the interaction is increased, eigenstates are strongly en¬ 
tangled and the local polarization is close to zero, its mi- 
crocanonical average. As predicted by ETH, each eigen¬ 
state is independently thermal and so the paradox of 
Quantum Thermalization is solved, as the memory of 
the initial condition fades out while entanglement grows 
through dephasing m- Instead, in the weakly interact¬ 
ing regime, for initial states close enough to exact eigen¬ 
states, a finite fraction of the polarization is doomed to 
survive [2TH23] . 

In absence of disorder, the spin-temperature is well- 
defined but very high m\ varying the ratio T//Acj e , the 
spin-temperature decreases up to a point where the sys¬ 


tem approaches the many-body localization (MBL), a dy¬ 
namical transition between an ETH and a non-ergodic 
phase Ell HO, surviving even in presence of microwaves 
[26] . In Fig. [4] (middle) we present a standard indicator 
for the transition: the variation of the local polarization 
between pairs of adjacent eigenstates versus the size of 
the system [4]. In the ETH phase, this quantity converges 
exponentially to zero indicating that all the fluctuations 
are more and more suppressed. On the contrary, in the 
localized phase, it saturates to a finite value, since fluctu¬ 
ations remain present even in the thermodynamic limit. 


Our results indicate that whenever the interaction with 
the environment is weak but not negligible, the dynamics 
reduces to quantum jumps between exact eigenstates of 
the electron system. Then, if ETH holds, the stationary 
state necessarily looks thermal, with few global param¬ 
eters (e.g. the spin-temperature) fixed relaxation and 
microwave irradiation. Instead, in the localized phase, 
only a weak DNP enhancement, triggered by few-body 
processes, can be observed. 


Spin temperature behavior. — It is important now to 
estimate the value of the spin temperature in the ETH 
phase. We first study a simplified model where the elec¬ 
trons in the EPR spectrum are assumed to be grouped 
into well separated macroscopic packets: 

N p 

H\,oy — ^ £> e + A p )J2S k z +vH in t (5) 

p k=l 

where ^2N p A p = 0 and N p is the number of electrons in 
the packet p. For 77 = 0 the spectrum of the Hamiltonian 
is composed by sectors of defined total magnetization and 
energy. The interactions are encoded in H int , which is 
chosen as a Gaussian random matrix inside each sector. 
When 77 is small - but still prevailing over the coupling 
with the bath - F int lifts the degeneracies in each sector 
selecting an ergodic basis in which the long-time density 
matrix is diagonal [27] . This model allows avoiding the 
numerical diagonalization, since statistical properties of 
the eigenstates are known. Moreover, the rates W n ^ n > 
in Eq. Q depend on 71, n' only via the matrix elements 
of local spin operators, | (n\ S l x \n') | 2 . Being the eigen¬ 
vectors perfectly ergodic in each sector, this quantity is 
actually determined by the pair of sectors containing re¬ 
spectively 71,77/, with weak statistical fluctuations m ■ 
These simplifications largely reduce the exponential dif¬ 
ficulty of the problem. In Fig [4] (right) we show that 
the stationary EPR spectrum for N = 64 spins perfectly 
agrees with the thermal Ansatz in Eq. ©• Moreover, /3 S 
and cjo in © can be fixed imposing that the energy and 
total magnetization become stationary for large times, 
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FIG. 4. Left: Density plot of the distribution of diagonal elements (n\ S l z \n) in the sector of vanishing total polarization. 
Colored regions represents, for each energy window (e, e + Ae), the smallest area containing half probability (U = 2 Mhz in 
blue, U — 15 Mhz in red, U — 45 Mhz in green). Middle: Logarithm of the variation of the local polarization between pairs 
of adjacent eigenstates of Eq.Q versus N . In the ergodic phase, this indicator vanishes exponentially in N . In the localized 
phase, it saturates to a finite value. Right: EPR spectrum for the toy model of Eq. © with 7 packets. The equilibrium profile 
is in blue. The yellow histograms show the stationary profile under MW irradiation with cjmw = 93.8775 Ghz and N = 64 
spins. The solid line is obtained from the Ansatz of Eq. © imposing the condition ©• 


which for the toy model leads to [19] 

2T 2 cJiP e (cjMw) + y ^N p ^ -- = 0 , (6a) 

p 1 

2T 2 ul A 0 P e (cj M w) + E -- = 0 . (6b) 

2Tl 

where uj p = uo e + A p , Po = — tanh(/3o; e /2) is the equi¬ 
librium polarization, Ao = uo e — wmw and we assumed 
that the microwaves only act on the resonating packet. 
Note that, for conserved quantities of ft as the energy 
and the total magnetization, the balance of the flows has 
a simple form since it reduces to the exchanges with the 
bath and microwaves. 

These results retrace the traditional prediction ob¬ 
tained within the phenomenological Ansatz of Eq. © 
proposed by Borghini [28]. Here, Eq. © naturally 
emerges once the strong suppression of fluctuations, char¬ 
acteristic of the ETH phase, has been assumed. However, 
the qualitative approach to hyperpolarization provided 
by this toy model largely underestimate the spin tem¬ 
perature value and hides its dependence from the micro¬ 
scopical parameters ([/, Ti,...) [29] . 

A richer scenario emerges instead from the exact diago- 
nalization of Eq. ([!]), where a stronger hyperpolarization 
enhancement is observed approaching the MBL transi¬ 
tion (see Fig. |3| and can be even amplified decreasing 
the relaxation time T\ (see Fig. [5|. Both effects agree 
with two well-known experiments, which fall beyond the 
applicability of the Borghini model [32]. The first showed 
that the enhancement occurs only at relatively low rad¬ 
ical concentrations PI, and therefore at weaker dipolar 
interactions. In the second, the addition of gadolinium 
complexes was used to induce a reduction of the relax¬ 
ation time T\ 133 HD- This, in turn, improved the signal 



FIG. 5. Ti shortening effect. The spin temperature /3 S , 
obtained from the fit of the electron polarizations (see Fig. [3j , 
is shown versus ljmw with U — 15 Mhz and different values of 
the relaxation time T \: the spin temperature is smaller when 
relaxation is faster, consistently with the observed increase in 
the hyperpolarization efficiency mm- 

enhancement and gadolinium is now commonly exploited 
in standard protocols for DNP sample preparation. 

Concluding remarks. — We presented a simple 
model for the study of DNP, providing a realistic depen¬ 
dence on the tunable parameters. The concept of out- 
of-equilibrium spin-temperature emerges naturally as a 
macroscopic manifestation of the ETH for the electron 
spin hamiltonian. 

Our study candidates DNP as a good ground for the 
direct observation of the MBL transition and its dynam¬ 
ical phase diagram. Two key advantages play in favor 
of this experimental setting. The first is that the two 
relevant control parameters for the transition are tun¬ 
able: U depends on the radical concentration and Auj e 
is proportional to the external magnetic field. The sec¬ 
ond is that the system does not require being isolated 
during the characteristic observation time, but, rather, 
that the relaxation is sufficiently slow to allow the pure 
quantum behavior to settle. Note that, in the past, 
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the spin-temperature polarization profile was already ex¬ 
perimentally observed in [33] for increasingly g-factors 
anisotropy, up to a critical value where the hole-burning 
profile popped out. The possibility of performing ex¬ 
periments precisely aimed at the observation of the elu¬ 
sive critical regime of the MBL is therefore concrete and 
promising. Moreover, the tunability of external parame¬ 
ters may allow the exploration of the phase diagram, even 
in regimes where the physics of spin glasses becomes rel¬ 
evant 34. 35]. 
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In the supporting information we provide the technical details of our computations. In Sec. I we write the master 
equation for the occupation probabilities of the exact eigenstates. In Sec. II we compute explicitly the polarization 
profile in absence of dipolar interactions, where P e (o;) displays an hole burning shape. In Sec. Ill, we derive the 
simplified master equation for the occupation probabilities of the sectors of the hamiltonian Htoy defined in the main 
text. Finally, in Sec. IV we specify the numerical procedures and the parameters used to produce the figures of the 
main text. 

The form of the transition rates, used in our calculations, is analogous to the one introduced by Hovav and 
collaborators, in Ref. [1]. In order to keep the paper self-contained, we derive, in the appendix, the relaxation scheme 
within the Lindblad formalism. In particular, the analysis of the energy exchange with the thermal bath is presented 
in App. A; the processes induced by the microwave irradiation are treated, within the rotating wave approximation, 
in App. B. Finally in App. C we show how to reduce the time evolution of the density matrix to a master equation 
for the occupation probabilities. 


I. MASTER EQUATION AND POLARIZATION PROFILE 

The spin dynamics is governed by a master equation for occupation probabilities, p n , of the 2 N exact eigenstates 
of H 


= E + w^] Pn , - [h( Unn ')wffi + w^f] Pn 

n'^n 


where the rate of microwave induced process writes 

/MW 4 CJ 2 T 2 


W ; = 
n ’ n l + (T 2 A ^ nn/ ) 2 


(?7/| S x |^)| 5 Aid n ,n' — ^n' (j> z ^z )^MW : 


(i) 


( 2 ) 


with e n , e n f and s™, s z ' eigenvalues of the operators H and S z on \n),\n') and S x = J2f=i $1 ( see App. B and App. C). 
In presence of a strong magnetic field: Acj 2 n , = (|e n — e n ' \ — u omw) 2 , as reported in the main text. The transitions 
with the thermal bath are governed by the function h(uj) = e fj +1 , which assures the detailed balance condition, and 
by the rates (see App. A) 


N 




(3) 


j = 1 

a=x,y,z 


The local polarization P e (uj, t ) of the electron spin under a magnetic field uo = uj e + A j writes as 

p e (uj, t)= 2 Tr[p(i)^] ~ 2 J^p n (t) (n\ S 3 Z \n). (4) 


where we approximated the density matrix p(t) with its diagonal. 


II. BLOCH EQUATIONS AND HOLE BURNING SHAPE 

In absence of dipolar interactions, the eigenstates of the spin system are products of local eigenstates of S J ZJ 

|a) = |tU---), (5) 

and the transition rates connect states that differ for a single spin flip, say the spin j : 

0J 2 iT 2 


W n 


_ w j _h(u e + Aj) 

^ + Aw 2 T 2 + 1 


(6) 



2 


where A uj = uj e + A j — cjmw- In this limit the probabilities p n factorize in the single spin probabilities, thus, the 
master equation (1) simplifies 

^= wi^ Pj=i (t) - w^ iPj=t (t) 

^ = W^ iPj=t (t) - w(^p j=l (t) 

and can be re-written as the Bloch equation for the local polarization P e (u:,t) = Pj=^(t) — Pj=i(t) at frequency 
(jj = u) e -f- Aj : 


dP e (u,t) _ 2 ujlT 2 P e {uj,t) P%(u) - P e (u,t) 

dt ~ Auj 2 T 2 +1 + Ti U 

where P®(uS) = — tanh(/?u;/2). The stationary solution corresponds to the hole burning shape of the irradiated EPR 
spectrum: 


(1 + A^ 2 T 2 )P>) 

1 + Ac u 2 T 2 + 2T 1 T 2 uj 2 


( 8 ) 


with Auj = uj — ujmw- 


III. MASTER EQUATION FOR THE TOY MODEL 


We now specialize the master equation (1) to the toy model defined by Eq. (4) of the main text. The spectrum of 
Ptoy? for V = 0, is highly degenerate and since [P t0 y: S z \ = 0, we label as V e}Sz the common eigenspace of Ptoy an d S z 
with eigenvalues, respectively, e and s z . When a small p > 0 is turned on, the matrix H[ nt resolves the degeneration 
and selects a non-degenerate basis in each subspace, so that for any | n) G V = V ejSz 

d 

i n ) = E c °i a ) ( 9 ) 

a= 1 

where the states | a) are product of local eigenstates of S \, as in (5). The sum runs over the states inside the subspace 
V, and d = dim(U). Choosing H int from the gaussian orthogonal ensemble (GOE), the coefficients c a are distributed 
according to the Porter-Thomas distribution [2]. In other words, for each eigenstate in a subspace V, they are chosen 
indipendently from the standard normal distribution, with the additional constraint of normalization 

E ( c «) 2 = 1 ( 10 ) 

4>=i 


so that c 2 = 1/d, where the average is over the random matrix ensemble. 

For p small enough, the energies of different sectors are well separated so that e remains a good quantum number. 
In this limit the master equation (1) connects only states belonging to subspaces differing by a single spin-flip. If we 
consider, for example, a spin i in the packet p that flips from up to down, we have a transition between a subspace 
V = V e ^ Sz and a subspace V' — V e -huj p ,s z - 1 , with uo p = uj e + A p . Taking | n) , \n f ) in V, V', the matrix element 
associated to this transition writes 

( n l &x \n') = c a c 'a ( a l SI |a'} = 1 E c «4 ( n ) 

aa' aGVj=t 


where we decomposed V = Vi—^ U Vi=i according to the value of the spin i and | a) is obtained from | a) by flipping i. 
From (11), we see that the matrix element converges to a gaussian number with zero average and variance 


I {n\&x W) I 2 


dim Vi-- 1 - 
4 dd' 


'Pjp ,t 
4 d! 


( 12 ) 
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FIG. 1. Sketch of the coarse-grained rate equation for the toy model defined in Eq. (4) of the main text. The rate of the 
transitions between two ergodic eigenstates \n) , \n') depends only on the two subspaces V, V'. The time-evolution of the 
system is then described by a master equation between subspaces. 


where d! = dim(V') and V p ^ is the probability that, for any state in V , a spin in the packet p is up. Neglecting 
the small fluctuations between states, we obtain that the rate of the transitions W n ^ n > depends only on the two 
subspaces that contain n and n', in this case W n ^ n > —> Wy^y. The master equation (1) between states can then 
be reduced to a coarse-grained dynamics between the occupation probabilities of each subspace, py — 5 ~2 n eyPn • The 
coarse-grained rates between subspaces write as (see Fig. 1) 


VIV-V' = d'Wy^y, = N p V p , t 


h(—hoj p ) 

Ti 1 + T 2 2 (cj p — u; M w) 2 


(13) 


Note that, the validity of Eq. (13) strongly relies on the ergodicity of the eigenstates. For instance, if the eigenstates 
were completely factorized, the rate of the transitions between | a) and | a') would retain finite fluctuations since the 
matrix elements | (a\ S l x \a') | 2 are either zero or one. 

Although the master equation described by Eq. (13) cannot be solved analytically, it allows the study of much larger 
systems as the number of subspaces grows only polynomially with N. Moreover, imposing that the d(H)/dt = 0 and 
d(S z )/dt = 0 , we obtain Eqs.(6) in the text. 


IV. METHODS 


The model defined by Eqs. (2,3) of the main text was solved with the following parameters 


Ti (sec.) 

T 2 (sec.) 

oji (27r Ghz) 

P (K- 1 ) 

U (Mhz) 

uj e (Ghz) 

Acj e (Mhz) 

0.25 -§- 4. 

HP 6 

0.25 x 10“ 4 

0.833 

2.0 4- 45.0 

93.9 

108.0 


and the local magnetic fields were chosen as 


Ai = Acj e 


2i - N - 1 
2 N 


(14) 


where N is the number of electrons. This choice corresponds in the limit of large N to a uniform distribution 
/(cj) ~ 1/Acj e . In order to obtain the polarization profile employed in Fig. 2, we performed the following steps: 

• Numerical diagonalization of H in each block at fixed total magnetization. 

• Solution of the linear system determining the stationary p n from the master equation (1). 
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• Evaluation of the stationary polarizations, P e ((o), using (4). 

The resulting polarization for each c j are then averaged over at least 50 realizations. 

For Fig. 3 (left), the distribution P of diagonal elements (n\ S\ \ n) was computed sampling, across i and 10 3 disorder 
realisations, all the eigenstates in the sector of vanishing total polarization. Then, the coloured regions were obtained, 
selecting for each energy window (e, e + Ae), the area with maximal concentration containing half probability. For 
Fig. 3 (middle), we focused on the sector of zero total magnetization, i.e. s z = 0, of H and computed the variation 
of the local polarization between pairs of adjacent eigenstates, defined as: 

= (n+l|S*|n+l)-(n|Sj|n) . (15) 

The expectation value \SM^\ was obtained by averaging over the spin index i , n belonging to the middle third of 
the spectrum, and a number of realizations ranging from 10 2 to 10 4 . For Fig. 3 (right), we determined numerically 
the stationary state of the coarse-grained master equation, defined by the rates in (13). We considered 7 packets 
resonating at frequencies uj p = cj e + Acj e (2p — 8)/14 with N p = 4, 8,12,16,12, 8,4. The probabilities V p ^ in (13) were 
computed by exact enumeration of the configurations of the subspaces. 


Appendix A: Lindblad equation in weak coupling approximation 

We summarize here the technique we use to treat the interaction of the spin system with the lattice bath. The full 
Hamiltonian is composed by three terms 

n = Hs + Hb + A^s-b (A1) 

where Hs is the spin system Hamiltonian, as given in Eq. (2) of the main text, Hb controls the evolution of the bath 
and Hs~b describes the coupling of the system with the bath; the parameter A measures the strength of this coupling. 
The most relevant contributions to Hs~b involve single-spin operators 

N 

Hs-b = 2 £ && (A2) 

3 = 1 

a=x,y,z 

where the operators <f>^ describe the excitations associated with the slow molecular modes localized around the 
electron spin j. We are interested in the effective long-time dynamics of the spin system when A is small, i.e. the 
interactions with the lattice are negligible with respect to the energy scales of the spin system. Within the Born- 
Markov approximation [3] the role of the bath is encoded in the autocorrelation function (4>^(to)4>^(to + £))bath> 
which being stationary depends only on t. For simplicity, we assume locality, homogeneity and isotropy so that we 
can restrict to a = /3, j = k and drop indexes. Then, we define the spectral density 

/ °° 

e iut Tr[pB®e- i 6 Bt $e i6Bt ] , (A3) 

-OO 

where pB is the density matrix of the bath. The quantum evolution for the density matrix of the spin system, p = ps> 
takes the form [4] 

d a 2 N 

^ = -i[H s ,p} + Cp=-i[H s ,p] + Y + ( A4 ) 

0 = 1 UJ 

a=x,y,z 

where C is the Lindbladian, which preserves positivity and trace of the density matrix. The sum over cj runs over 
all the energy differences of the exact eigenstates of Hs. The presence of random dipolar coupling in Hs ensures 
non-degenerate gaps so that 


i^nn' ) — \n ) (jl \ S ^ \n) (jl\ , CO nn ' — 6 n C n ', (A 5 ) 

with e n , e n r the energy levels associated to the eigenstates \n),\n'). Let us observe that in the weak coupling limit, the 
bath influences the dynamics by allowing transitions between the exact eigenstates of the system Hamiltonian Hs- 
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Since the bath is at equilibrium, it is easy to check, by replacing t t — i/3 in (A3), the KMS condition [5, 6] 


J(lj) = e^J(-co) . (A6) 

This condition implies the detailed balance of the rates and therefore the convergence toward equilibrium p = e~^ Hs jZ 
where yd -1 is the bath temperature. 

It is convenient to rewrite Eq. (A4) separating the evolution of diagonal and off-diagonal terms of p. For the 
diagonal terms, we have 


dpnn 

dt 


^ ^ [h( W nn ') Pn'n' hirin')Pnn] W nn ' 

n' 


(A7) 


where h(ui) 


e /3u, + l 5 


satisfying (A6), and the rate of the transitions between the levels n and n' is given by 


W hat J 1 

v v n,n' 


2 

Ti 


N 

E K n l^al n ')| 2 - 

3 = 1 

a=x,y,z 


(A8) 


The relaxation time, Ti, fixes the time scale of the electronic single-spin transitions. Its value depends on J(uj) and 
the coupling strength A. The transitions involving S XiV correspond to a change of the total polarization and their 
time scale can be estimated by the experimental value of the electron relaxation time, T\ ~ 1 sec. On the contrary 
a direct measure of the transitions involving S z is difficult and they are usually neglected [1]. Here we assume that 
their characteristic time scale remains of the order of 1 sec. For the off diagonal terms we have 


dp n 


dt 


— i (e n Cm) Pnrn 


Pnrri 

2~2,nm 


(A9) 


The exact expression of T 2 ?nm i n terms of the spectral function can be derived from Eq. (A4). However since the off 
diagonal decay rate is generically much faster then the relaxation time T\ we assume for simplicity T 2 ?nm = T 2 , with 
the experimental decoherence time T 2 = 10 -6 sec. [1]. 


Appendix B: Microwave irradiation in rotating wave approximation 

We irradiate the system with a microwave field 

^mw = 2cji cos(cjmw t ) S x (Bl) 

where uj\ is the intensity of the field, ujmw its frequency of the order of 100 Ghz and S x = S 3 X . The evolution 

equation (A4) is modified as 

^ = -i[H s ,p}-i[H M w,p}+Cp. (B2) 

Since the microwave power is much smaller than the energy scales of the spin system, the Lindbladian, £, is not 
affected by the microwave field. Eq. (B2) is now time-dependent and it is convenient to employ the so-called rotating 
wave approximation (RWA). We define the density matrix in the rotating frame as p^ = U(t) p£H(t), with U(t) = 
exp(iS' 2 co’MW^) and S z = JT &Z- ^^is way, since [Hs,S z ] = 0, Eq. (B2) writes 

^ - i[U(t)H M wUHt),P {r) } + U(t ) (Cp) U\t) (B3) 

where = H$ — wmw S z . Moreover we have 

N 

^MW = U(t)H MW U\t) = ~ wiS* 

3 = 1 
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where the fast oscillating term has been neglected. Since the total magnetizanion S z is a good quantum number, we 
have from (A5) 


U(t)Si(±u)U\t) = e ±i ( s "- s ” ,)u ' Mwt S , ^(±w) 


(B4) 


where s”,s” are the eigenvalues of S z on \n),\n'). It follows that the last term of Eq.(B3) can be rewritten as £p tr> . 
and in conclusion, adopting the RWA approximation, Eq. (B2) becomes 


^ -i^[s x ,p {r) ] + £p (r) ■ 

Separating the diagonal and off-diagonal dynamics, it can be rewritten as 


7 (r) 

dpnn 

dt 

J r ) 


= E H Sx ( Pn'n ~ Pnn') + E h (~ u nn’ K'n' ~ h(u} nn >)p^ 


Tx/'bath 

W n,n' 


= - fiAuj nm + -A - iU! y] {n\ S x \ ri) - (n'| S x I to) pp n , 


dt \ T 2 

where Aui nm = e n - e m - (s” - s™)w M w- 


(B5) 


(B6) 


Appendix C: Master equation in the fast dephasing limit 

In our system, the dephasing time T 2 is much faster than the relaxation time T\ (T 2 /T 1 10 -6 ) and the quantum 

evolution can be reduced to a master equation for the diagonal elements of the density matrix in the basis which 
diagonalises Hs • We start rewriting (B6) as 


^ r =L 0 p (r) +L 1 p^) (Cl) 

where Lo,i are superoperators acting linearly on the space of densitiy matrices. Lo is the dominant part involving 
contributions from T 2 -1 and Aco> nn /, while Li contains the remaining contributions. We denote as e nn / the matrices 
with a single one in position n, r\!. They correspond to the right eigenvectors of Lq with eigenvalues: 


^O^-nn' 


( ~(T 2 1 + iAu nn r) n^n ’, 
| 0 n = n '. 


(C2) 


In the long-time limit the dynamics is restricted to the subspace with eigenvalue 0, corresponding to the diagonal 
density matrices. The effective evolution in this subspace due to the presence of Li can be derived by standard 
perturbation theory. By setting p n = pnn = p nn and indicating as (Li) nn q mm / the components of the superoperator 
Li on the basis e nr7 /, we obtain 


dp n 

dt 


E 


(Li), 


E 

rn^rn' 


(Ll)r 


'(Li) 


mm' ,n'n' 


^ 2 , mm' + iAcUn 


Pn' 


After some algebra we get the final expression given in Eq.(l) 


(C3) 
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